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Abstract: We study F(R) modified gravity models which are capable of driving the 
accelerating epoch of the Universe at the present time whilst not destroying the standard 
Big Bang and inflationary cosmology. Recent studies have shown that a weak curvature 
singularity with — )• oo can arise generically in viable F{R) models of present dark energy 
(DE) signaling an internal incompleteness of these models. In this work we study how this 
problem is cured by adding a quadratic correction with a sufficiently small coefficient to the 
F(R) function at large curvatures. At the same time, this correction eliminates two more 
serious problems of previously constructed viable F(R) DE models: unboundedness of the 
mass of a scalar particle (scalaron) arising in F{R) gravity and the scalaron overabundance 
problem. Such carefully constructed models can also yield both an early time inflationary 
epoch and a late time de Sitter phase with vastly different values of R. The reheating epoch 
in these combined models of primordial and present dark energy is completely different from 
that of the old R + R'^/GM'^ inflationary model, mainly due to the fact that values of the 
effective gravitational constant at low and intermediate curvatures are different for positive 
and negative R. This changes the number of e-folds during the observable part of inflation 
that results in a different value of the primordial power spectrum index. 
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1. Introduction 

Numerous recent observational data prove convincingly that the Universe is undergoing 
accelerated expansion at the present time, whilst decelerating in the past for redshifts larger 
than about z ~ 0.7. If interpreted in terms of the Einstein general theory of relativity, 
this acceleration requires the existence of some new component in the right-hand side of 
the Einstein equations, dubbed dark energy (DE), which remains practically non-clustered 
at all scales at which gravitational clustering of baryonic and dark non-baryonic matter is 
observed, and which has an effective pressure pz)^; approximately equal to minus its effective 
energy density poE- Thus, its properties are very close to those of a cosmological constant 
A (see [1, 2, 3, 4, 5] for some reviews). The simplest possible DE model, A combined with a 
non-relativistic non-baryonic dark matter (the standard spatially flat ACDM cosmological 
model), is completely self-consistent from the mathematical point of view and provides 
a good fit to all existing observational data [6]. In this case A acquires the status of a 
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new fundamental physical constant. However, its required value is very small compared to 
known atomic and elementary particle scales (which are well below the Planck scale), so a 
firm theoretical prediction for this quantity from first principles is currently lacking. 

On the other hand, in the second case when a component with qualitatively similar 
properties is assumed to exist - in the inflationary scenario of the early Universe, we are 
sure that this "primordial DE" may not be an exact cosmological constant since it should 
decay in the early Universe. Hence it is natural to seek non-stationary models of the current 
DE, too. 

Among them, the simplest purely gravitational models in 3 + 1 space-time dimensions 
are provided by F{R) gravity which modifies and generalizes Einstein gravity by incor- 
porating a new phenomenological function of the Ricci scalar R. F{R). They represent 

a self-consistent and non-trivial alternative to the ACDM model. The literature on these 
models is dense, and we direct the reader to [7] and references therein for a detailed recent 
review. The action of F(R) gravity is given by 

S = ^ y d^x^FiR) + Sm , (1.1) 

where Mpj = l/SvrG, ^ = c = 1 is assumed throughout the paper, and Sm describes all 
non-gravitational matter including non-relativistic (cold) dark matter which is minimally 
coupled to gravity.^ The field equations following from (1.1) have the form 

F'{R)Rl - \f{R)5^^ + {5^^U - V,VnF'{R) = Mpf t; , (1.2) 
and their trace reads 

3DF'{R) + RF'{R) - 2F{R) = M'^T , (1.3) 

where T is the trace of the matter energy-momentum tensor and the prime denotes the 
ordinary derivative with respect to an argument. 

In the Jordan frame where the action (1.1) is written, fermion masses are constant 
and atomic clocks built from usual matter measure the proper time. However, equivalent 
description is possible in the Einstein frame where gravity resembles standard GR but free 
particles of usual matter do not follow space-time geodesies due to an interaction with a 
new scalar field. As will be seen below, sometimes it is easier to solve equations in the 
latter frame first. We also assume the metric variation of (1.1); the Palatini variation of 
formally the same action leads to a completely different theory, in which the number of 
degrees of freedom is not the same. 

In the absence of matter, exact de Sitter (positive constant curvature) solutions of Eqs. 
(1.2) are given by real positive roots of the functional equation 

RF'{R) - 2F{R) = . (1.4) 

^The sign conventions are: the metric signature (-- 1— I— 1-), the curvature tensor R'^^pi, = dpV^^ — ..., Rf^^ = 
R%u, so that the Ricci scalar R = R'j^ > for the de Sitter space-time and the matter-dominated cosmo- 
logical epoch. 
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These solutions (and solutions close to them) are the basis for a description of primordial 
and present DE. They are future stable if 



F'{Ri)/F"{Ri) > Ri 



(1.5) 



where Ri is a root of Eq. (1.4). This condition was first obtained in [9], the easiest way to 
derive it is to use the trace equation (1.3) for a small perturbation R— Ri (in fact, F"{Ri) 
should be positive, too, as will be discussed below). 

F{R) gravity is a special class of scalar-tensor gravity with a vanishing Brans-Dickc 
parameter ojbd- If F"{R) is not zero identically, it contains a new scalar degree of freedom 
dubbed "scalaron" in [8], thus, it is a non-perturbative generalization of Einstein gravity. 
We will consider it as a purely phenomenological semiclassical macroscopic theory of gravity 
which arises from some more fundamental quantum microscopic theory after tracing out 
degrees of freedom which are not excited at sufficiently small space-time curvature. Thus, 
the resulting function F{R) need not necessarily be some simple (e.g. polynomial) function 
of R. It may well have some complicated behaviour for small R, too, as many examples 
from condensed matter physics teach us. So, we will not discuss which functional form of 
F{R) is "natural" in any sense. On the other hand, for a phenomenological F{R) model 
to be viable, it should satisfy a rather large list of viability conditions: 

1) Classical and quantum stability in the region of R where we want to use this theory: 



The first condition means that gravity is attractive and the graviton is not a ghost. It was 
recognized long ago that its violation during the time evolution of a Friedmann-Robertson- 
Walker (FRW) background results in the immediate loss of homogeneity and isotropy and 
formation of a strong space-like anisotropic curvature singularity [10, 11]. The second 
condition on the flat background was also known since the first papers on F{R) gravity 
[12], and was assumed when constructing inflationary models in F{R) gravity [8]. However, 
in the case of F{R) models of present DE, the necessity to keep it valid for all values of R 
during the matter- and radiation-dominated stages in order to avoid the Dolgov-Kawasaki 
instability [13] has been realized rather recently [14, 15]. In addition, a weak ("sudden") 
curvature singularity forms generically if F"{R) becomes zero for a finite value R = Rg. 
This is also undesirable; see the discussion below in Sec. 2.4. 

2) Existence of the stable Newtonian limit for all values of R where Newtonian gravity 
accurately describes observed inhomogeneities and compact objects in the Universe, i.e. 
for R exceeding the present Priedmann-Robertson- Walker (FRW) background value Rq = 
R{to), where to is the present moment, and up to curvatures in the centre of neutron stars. 
The conditions required for this are 



for R ^ Rq. Note that in this regime the effective scalaron mass squared is AIg{R) = 
1/{3F"{R)), as directly follows from Eq. (1.3). Then the second of the conditions (1.6) 
means that the scalaron is not a tachyon, while the last of the conditions (1.7) implies 



F'{R) > 0, F"{R) > . 



(1.6) 



\F{R) -R\<.R, \F'{R) - 1| < 1 , RF"iR) < 1 



(1.7) 
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that its Compton wavelength is much less than the radius of curvature of the background 
space-time. For more general backgrounds than matter-dominated FRW, for which General 
Relativity (GR) has to be used in full (in particular, if the pressure P of matter is not small 
compared to its energy density p), the conditions (1.7) guarantee that non-GR corrections 
to a space-time metric remain small. 

3) Absence of deviations from GR at the level of accuracy following from present 
laboratory and Solar system tests of gravity. 

4) Existence of a future stable (or at least metastable) de Sitter asymptote. This is 
necessary for a description of the present DE, which behaves in a similar manner to that 
of a cosmological constant. 

5) F(R) cosmology should not destroy previous successes of present and early Universe 
cosmology in the scope of GR including the existence of the matter-dominated stage driven 
by non-relativistic matter preceded by the radiation-dominated stage with the correct Big 
Bang nucleosynthesis (BBN) of light elements and, as we shall see, some kind of inflation 
prior to this. 

The first of these conditions is also the main reason why we do not include other 
invariants constructed from the Riemann tensor and its derivatives as additional argu- 
ments of the function F. Indeed, it has long been known [16, 17] that if F also depends 
on R^yapR^'^'^^ then gcncrically a new massive spin-2 particle appears which is a ghost 
whenever the standard masslcss graviton is not a ghost. This is problematic from the 
quantum field theory point of view for many reasons, see e.g. [18]. Moreover, as argued in 
[19], cosmological models with ghosts are unsatisfactory even at the purely classical level. 
In particular, in such models we can no longer explain the observed approximate large- 
scale homogeneity and isotropy of the Universe without tremendous fine-tuning of initial 
conditions, even if we include a primordial inflationary stage. The only way to avoid this 
ghost (without including derivatives of the Riemann tensor in the action) is to consider 
F = F{R,G) where G is the Gauss-Bonnet invariant, G = R^uapR'^''"'^ - ^RnuR'^'' + R^ ■ 
However, it was recently shown in [20] that in this case linear scalar perturbations have 
pathological behaviour in the ultra-violet regime, except in some special cases. ^ Invariants 
in F containing derivatives of the Riemann tensor lead to new particles (in particular, scalar 
ones if F depends on R and its derivatives only [22]) among which ghosts are generic, too, 
see e.g. [23]. Thus, such terms may be considered when making perturbative expansion 
around solutions of the Einstein gravity, but (possibly apart from some exceptional cases 
still to be found) they are of no use in our approach since we want to use a (semi) classical 
modified theory of gravity in a fully non-perturbative regime. 

For these reasons, we concentrate on F{R) effective macroscopic gravity only. Note 
that this is in contrast to quantum-gravitational and string corrections to Einstein gravity, 
which generically produce terms with all possible invariants of the Riemann tensor. The 
required "i?-dominance" presents a serious problem for such microscopic mechanisms to act 
as the origin of F{R) gravity, however there exists a number of cases when just this form of 

^After the first variant of our paper was submitted to JCAP and archives, the further paper [21] appeared 
where this remaining special case was shown to possess an ultra-violet instabiUty of scalar perturbations 
on a FRW background with matter in the form of a perfect fluid, too. 
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modified gravity appears in some limit. The simplest of them follows from above-mentioned 
fact that F(R) gravity is a particular case of more general scalar-tensor gravity with the 
Brans-Dicke parameter ujbd = 0. Therefore, it also yields a good approximate description 
for scalar-tensor gravity with I^bdI ^ 1, in particular, for a non-minimally coupled scalar 
field with a negative and large by modulus coupling constant ^ (we use the sign convention 
where conformal coupling corresponds to = +1/6), cf. [24]. For this reason, in particular, 
predictions of the Higgs inflationary model [25] (without loop corrections to the Higgs 
potential) for primordial power spectra of scalar (density) perturbations and gravitational 
waves are the same as those of the F{R) = R + i?^/6M^ model - the simplified variant of 
the model [8]. Another, completely unrelated case where phenomenological F{R) gravity 
arises rather unexpectedly [26] is the so called emergent gravity approach using ideas and 
methods borrowed from quantum theory of condensed matter. This illustrates the well- 
known fact that an elegant and internally consistent mathematical model may appear 
multiple times from totally different physical foundations. Note that in both of these 
examples, the new scalar gravitational degree of freedom (scalaron) is present already at 
the underlying microscopic level; however, this does not mean that it is fundamental even 
at this level. 

Although F{R) gravity can successfully pass the first requirement from the list above, it 
is evident from the beginning that only a very narrow subset of all possible F{R) functions 
may be of interest for cosmology. The situation is simpler in case of inflationary models. 
Here, the simplest variant with F"{R) ^ identically, namely F = R + R? /QM'^ (where M 
is the scalaron rest-mass at low curvature), presents an internally self-consistent inflationary 
model with slow-roll during inflation (see Sec. 4.4 below) and a graceful exit to a subsequent 
FRW matter-dominated stage driven by scalarons [8] . Reheating, creation of usual matter 
and transition to the radiation-dominated FRW stage are achieved by gravitational particle 
production due to strong oscillations of R during the scalaron-dominated stage [8], see 
[27, 28, 29, 30, 31] for more details regarding the background evolution and reheating in 
this model. In addition, in contrast to many other inflationary models proposed later which 
have been falsified by observational data, this model still remains viable since it predicts 
the value of the slope of the primordial power spectrum of scalar (density) perturbations 
ng = 1 — 2/N [32, 33] in agreement with present observational data. Here A'' is the number 
of e-folds between the first Hubble radius crossing of the present inverse comoving scale 
0.002 Mpc~^ and the end of inflation; N ^ (50—55) for the reheating mechanism mentioned 
above. The tensor-to-scalar ratio is rather small, r = 12/Ar^, but not negligible [33]; see 
also [34, 35, 36] for further papers on equations and solutions for perturbations and [37] for 
the energy theorem in this model. To fit the observed amplitude of the power spectrum, 
the only free model parameter M should be chosen as M ~ 1.5 x 10~^(A/50)^^Mpi. 
Furthermore, we will show in Sec. 4.4. that all viable inflationary models in F{R) gravity 
with other values of Ug and r have behaviour close to R^ at large R (or around some large 
fixed value R = Ri). Thus, the R^ behaviour of F{R) is characteristic for inflation in F{R) 
gravity. 

Due to the remarkable quantitative analogy between properties of primordial DE sup- 
porting inflation in the early Universe and present DE, it is tempting to use F(R) gravity 
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to also build models of the latter, to act as alternatives to the trivial case of a cosmological 
constant which corresponds to F{R) = R — 2A. Indeed, many such models have been 
proposed, beginning with the papers [38, 39, 40, 41]. However, most of these attempts 
remained unsuccessful since the conditions written above had not been fully satisfied (once 
more, we direct a reader to the review [7] for an extensive list of papers on the subject). 
As a result, there has arisen a widespread doubt if viable F[R) models of present DE may 
exist at all (in contrast to primordial DE). 

Still at last a rather narrow class of functional forms of F[R) was found [42, 43, 44] 
which can satisfy the first four viability conditions from the list above and even partly the 
fifth one, with regards to the existence of a stable matter dominated epoch with a{t) oc t^/^ 
in the recent past which is driven by cold dark matter and baryons (not, we stress, by terms 
due to modified gravity.) Here, a{t) is the scale factor in a FRW spacetime. It should be 
emphasized from the very beginning that in none of these viable F{R) DE models is it 
possible to derive the energy scale of present DE from first principles. Instead, it has to be 
inserted into the action (1.1) as a free parameter, the value of which is taken from observa- 
tional data. Thus, F{R) DE models may not be superior to DE being an exact cosmological 
constant A, they are simply an alternative to it. In the models [42, 43, 44], F{R) is analytic 
for R = 0, has a non-trivial structure for i? ~ i?(to) and then quickly approaches Einstein 
gravity with an effective cosmological constant for R ^ Rq, sec Eq. (2.3) below. Also, in 
these models the condition F{0) = (dubbed the "disappearing cosmological constant" 
in [44]) is imposed by hand to ensure that there is no true cosmological constant in flat 
space-time (otherwise, why work with F{R) gravity at all?) The main difference between 
these models is in the law of approach to the standard gravity for i? — >■ oo: an inverse 
power law in [42, 44] and an exponential in [43]. 

However, it was immediately recognized that the story of constructing at least one 
viable cosmological model of present DE in F{R) gravity was not finished: three new 
problems related to the last, fifth viability condition arise when tracing small deviations 
of a FRW background from the standard ACDM model to the past. First, as was shown 
in [44] (see also [45]), for these models the frequency of small oscillations of the Ricci 
scalar R (i.e. the scalaron rest-mass) around the general relativistic limit iZcR = ~F/Mp^, 
Lu = Ms{R) = 1 / "iF" {R) , grows quickly to the past, i ^ 0, and exceeds the Planck value 
very soon, thus invalidating classical consideration of the theory (1.1). For example, for the 
model [43] it happens already for redshifts z > 1 and matter densities exceeding ~ 10~^^ 
g cm~^. 

The second problem is that the amplitude of these linear oscillations quickly grows back 
in time [44]. As a result, we get the scalaron overabundance problem which is actually 
a problem of initial conditions of the Universe; they should be such that the scalaron 
number density (which is proportional to the square of an amplitude of these oscillations) 
should be sufficiently small at the period of BBN. Third, due to the same reason, linear 
consideration of these oscillations may become inadequate even before their backrcaction 
on a FRW background becomes important. A non-linear analysis of the these models was 
undertaken in [46], where it was observed that the Ricci scalar would gencrically evolve to 
a weak singularity at some finite time in the past. This singular behaviour was predicted 
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independently in [47] for a more general class of models. These three problems once again 
raised the question as to whether viable DE F(R) gravity models can be constructed. 

The main aim of this paper is to find a way to avoid these three difficulties and to show 
that there indeed exist F{R) DE models satisfying all five viability conditions. We shall 
show that for this purpose, as envisaged in [44], it is sufficient to change the behaviour of 
F{R) at large R by adding the term oc R^ with a sufficiently small coefficient to ensure 
the existence of an inflationary stage in the early Universe. Moreover, since this behaviour 
is characteristic for viable inflationary models in F(R) gravity, we show additionally that 
it is possible to construct a combined model where both primordial and present DE are 
described in the scope of F{R) gravity.^ Rather unexpectedly, such a construction also 
requires us to change the low curvature behaviour of F{R) in the models [42, 43, 44] for 
R < Rq and further to the region of negative values of R which are not observable at the 
present time, in order to avoid violation of the conditions (1.6) during strong oscillations 
of R after the end of inflation. Moreover, it will be shown that a non-trivial structure 
of F{R) at low R, required for an alternative description of present DE (as opposed to a 
cosmological constant), greatly affects the stage of post-inflationary evolution and reheating 
in this combined model, and even results in the change of numerical values for parameters 
of primordial power spectra of perturbations generated after inflation. 

The rest of the paper proceeds as follows. In section 2 wc review the cosmological 
evolution of viable models of present DE in F(R) gravity, and comment on the existence and 
properties of the weak singularity mentioned above. We then consider possible approaches 
to eliminate this singularity and bound the scalaron mass by introducing additional terms 
into the action in section 3. Finally, we consider the possibility that these additional 
terms may drive an early inflationary period of the Universe, and study the slow-roll and 
reheating epochs in section 4. Section 5 contains conclusions and discussion. 

2. Review of cosmological evolution 

We will be concerned with the following F{R) functions describing present DE which satisfy 
the first four viability conditions listed above and possess a stable matter-dominated epoch 
at intermediate redshifts z for some range of their parameters: 



where b and c are dimensionless constants, eAB = Rvac/{b + log(2cosh6)) and n > 0. The 
model (2.1) is the model introduced in [42] (with a slightly different notation of parameters); 

^Of course, this model does not really unify primordial and present DE since we have to introduce two 
tremendously different curvature scales corresponding to inflation in the early Universe and to the present 
space-time curvature by hand. 




(2.1) 



(2.2) 
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it is similar to the model of [44]. The model (2.2) is from the paper [43] (a model with 
a similar behaviour was also introduced later in [45]). For large R, these models mimic 
General Relativity with a cosmological constant, in the sense that for R S> -Rvac, -^HSS and 
Fab can be expanded in the region of interest to cosmology as 

FiR)^R-^ + xiR) . (2.3) 

Thus, -Rvac/4 acts as a small effective cosmological constant induced by space-time cur- 
vature. Since F(0) = 0, there is no true cosmological constant in these models. x(i?), 
x'iR) and x"iR) arc all small functions of R, which satisfy x(ii)/i2 <C 1, x'(-R) ^ 1 and 
Rx"{R) < 1 for i? > i?vac- For the HSS and AB models, we have 



2n+l 

XHSS = XAB = -^e^' exp (-2i?/eAB) , (2.4) 



where eAB and enss = ^vac/(2c)^/^^"+^' arc smaller than i?vac- 

In this section we briefly review the cosmological evolution of the AB and HSS models, 
in particular the behaviour of the scalar degree of freedom (scalaron), and highlight the 
existence of a singularity in the evolution of the Ricci scalar. In ref.[44], the trace equation 
(1.3) was solved using a perturbative approach: an ansatz R = RQji + 6R was taken, where 
Rqji = —T/Mp^, and the equation linearized for SR <C Rgr- It was found that the Ricci 
scalar for the HSS model will generically undergo rapid oscillations around its General 
Relativistic limit, a result also obtained for the AB model in ref.[45] (although see ref.[46] 
for a discussion of the applicability of the linearized approach in the case of the AB model.) 
The frequency and amplitude of these oscillations increase without bound for both models 
as R grows to the past (t — > 0). Specifically, SR contains an oscillating component that is 
given by 



SRosc = Ca-3/2(F"(i?GR))-3/4sin 



dt 



V3F"(i?GR) 



(2.5) 



It represents scalaron oscillations (particles) with the frequency (rest-mass) Ms(R) = 
(3F"(/Zgr))~^^^ in the regime when > \R\gr- Typically, Mg becomes » Mpi during 
the matter and radiation eras for both models. 



2.1 Non-linear oscillations and existence of a "sudden" singularity with — )• oo 

In refs.[46], an alternative approach to solving (1.3) was considered. In this work it was 
noted that at curvatures of cosmological interest, we can use the expansions (2.4) to write 
(1.3) as a non-linear oscillator equation. By making the field redefinitions 



R = Rgk- ^iog{l + x), 



R 



Rgr 



(l + a;)V(2n+i)' 



(2.6) 
(2.7) 
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for the AB and HSS models respectively, and neglecting terms of order 0{x{Rgr)), (1-3) 
can be written as 



XAA + ^XA+|log(l + i)«0 (AB), (2.8) 

x« + ^XA-|(^j^-l).0 (HSS), (2.9) 

for R ^ i?vac, where we have taken n = 1 for simplicity in the HSS model. In (2.8) and 
(2.9) a; is a dimensionless field, e has the same dimensions as the Ricci scalar, and a{X) and 
P{X) are dimensionless damping terms that will be unimportant in the following discussion. 
A subscripts denote derivatives with respect to the 'fast time' A, which is related to the 
cosmological time by 

' (2.10) 



where we have introduced the dimensionless function /q = F'(i?GR) — 1- Both (2.8) and 
(2.9) are damped, non-linear oscillator equations, with potentials given by integrating the 
expressions dV/dx = (e/6) log(l + x) and dV/dx = (e/6)(l - (1 + x)-V3) for the AB and 
HSS models respectively. 

Equations (2.8) and (2.9) can be solved numerically, and R obtained from (2.6) and 
(2.7). By taking into account non-linear terms in (1.3), a number of effects were observed 
that were not apparent in the linearized analysis. Specifically, it was found the oscillations 
become asymmetric in the past (as expected, when the linearized approximation breaks 
down), and further that by evolving R backwards in time, it will generically evolve to a 
singularity. 

A singularity arises in both models as a: — t- — 1 after a finite time. At this point, 
the potential has the asymptotic behaviour V{x) — )• const, but its derivative diverges 
dV/dx oo. Prom (2.6) and (2.7), it is clear that the Ricci scalar will diverge as x ^ — 1. 
The potentials for the AB and HSS models are given by 



VMi{x) = ai + -(1 + x) [log(l + x) - 1] , 



^Hss(a;) = "2 + g 



(l+^)_^!L±l(l + a;)2n/(2n+l) 



(2.11) 
(2.12) 



and are exhibited in fig.l. ai,2 are integration constants, which dictate the ground state 
of X. In the limit x — > — 1, we have V const, as expected. 

The divergence of the Ricci scalar represents a weak singularity in the sense that a 
and R diverge whilst p, p and a remain finite. As a result, a small vicinity of space- 
time containing this space-like singularity (and from both sides of it) may be covered by 
the Minkowski metric with small perturbations which, however, are not continuous. 
Moreover, if the singularity occurs during the matter-dominated epoch when p <^ p, the 
Newtonian approximation is applicable around it. Let us obtain the expression for the 



-9- 



V(x)/e 




1+x 



1 2 3 4 5 6 

Figure 1: The potential (2.11) for the AB model (top), with 6 = 4, Rgr = 10i?vac and ai = 0, 
and HSS model (2.12) (bottom), with n = 2, c = 1, i?GR = 10i?vac and a2 — 0. The singular point 
is at a; — >■ — 1, where the potential is finite but dV/dx diverges. 

behaviour of a FRW scale factor near the singularity. It follows from the trace equation 
Eq. (1.3) that the analytic form of a{t) in the vicinity t = tg - the time at which R 
diverges - is given by 

a{t) = ao + ai{t- ts) + a2{t - t^f (log |t - ts| + Sa) + ... , (2.13) 
a{t) = ao + ai(t-t,)+a2|t-ts|^'+'"^/('+") + ... , (2.14) 

for the AB and HSS models respectively, where ao,i,2 and 0,2 are constants. Such a sin- 
gularity was considered in [48] from a kinematic viewpoint (i.e. not as a solution of any 
dynamical equations) where it was called "sudden" . It also appeared in a different dynam- 
ical setting in [49], where it was dubbed the "Big Boost". 

In contrast to strong curvature singularities occurring in cosmology and inside black 
holes, there is no geodesic incompleteness here [50]. However, this does not mean that the 
weak "sudden" singularities (2.13, 2.14) are harmless. Just the opposite, they are unde- 
sirable. Indeed, it can be checked that Eqs. (1.2) do not supply us with any information 
on how the coefficients 02 for t < tg and t > tg are related to each other. Thus, the "sud- 
den" singularity results in the loss of predictability; from initial data given at any Cauchy 
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hypersurface with t = const < ts, it is not possible to predict the space-time metric for 
t > ts unambiguously. Note also that in contrast to strong curvature singularities, we may 
not evade the problem by arguing that quantum-gravitational effects invalidate the very 
notion of deterministic classical space-time near such a singularity. It can be shown [51] 
that such effects remain small as one approaches the weak singularity; the energy density 
and other components of an averaged energy-momentum tensor (EMT) of quantum fields 
in such curved space-time generically remain subdominant to the background (this also 
follows from general expressions for an EMT average value in a weakly curved background 
obtained in [52]). So, the singularity has to be resolved at the classical level. 

Therefore, the appearance of "sudden" singularities in a given F{R) DE model signals 
the internal incompleteness of said model. One way to resolve the singularity is to return 
from the effective macroscopic F{R) theory (1.1) to the underlying microscopic theory from 
which the former originated, and sec what happens in the latter. Then, however, there is 
no guarantee that the resulting ansatz for the resolution of the "sudden" singularity will 
be universal and will not depend on the choice of underlying physics. Since we prefer to 
remain at the purely phenomenological level in this paper, we instead look for a way to 
modify the given F{R) model in such a way that the weak singularities do not appear 
at all, at least in solutions which are of interest for cosmology, an approach that we will 
consider in Sec. 3. 

2.2 General viable F{R) models of present DE 

Thus far, we have focussed attention on two specific models. We will now generalize 
the existence of a singularity to an arbitrary F{R) function for which F"{R) > and 
R\F'{R) — 1| — as i? — )■ oo. Whilst viable F{R) models do not have to satisfy this 
asymptotic limit at infinity, they do have to satisfy RF"[R) <^ 1 for R S> Rva.c, as we 
demand that any modification to General Relativity should only become dominant at late 
times. We begin with the function F{R) = R — -Rvac/2 + x{R): where x{R) satisfies 
X{R)/R < 1, x'iR) < 1 and Rx"{R) < 1- We then define 1 + x = x'{R)/Xo, where 
Xq = x'(-Rgr)) and write (1.3) in terms of x, 

3nxo(l + x)-{R- Rgr) + XqR{^ + 2x{R) ^ 0, (2.15) 

where we have defined Rqr = —T /M'^^+R^g^. It should be understood that R = R{x, -Rgr) 
in (2.15), since the Ricci scalar will generically be a function of both x and -Rqr. Using a 
flat FRW metric ansatz, we obtain 



'r _L Ton '^ _LQH-1 l_Ln_L M-i / , R{x,Rgr)) - Rgr , '^xjR) - R{x)Xoi'^ + x) 
Xo [xxx + [2(logXo)A + 3-ff] xx\+{l+x)nxo+ h = 0, 

(2.16) 

where H = ax/a. Now by using the conditions x'{R) ^ 1 and x{R)/R ^ 1) Eq. (2.16) 
can be approximated as 

X'o [xxx + [2(log x'o)x + 3^] xx] + ^ 0. (2.17) 
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The x'o function multiplying the xaa and x\ terms is unimportant and can be removed by 
a redefinition of the time coordinate. The final step is to invert the expression 1 + x = 
x'{R)/Xo to obtain R = R{x, Rqr) and substitute this in (2.17). Then we can associate the 
last term on the left-hand side of (2.17) with a potential gradient dV/dx = {R{x,Rgr) — 
i?GR)/3. 

The two models considered in this paper have particularly simple xi^) functions, 
and hence it is straightforward to invert the expression 1 + x = x'(-R)/Xo fo obtain R = 
R{x, Rqr), giving an analytic form for dV/dx. However, in general we will not be able 
to write dV/dx in terms of known functions. Another feature of the two models under 
consideration is that all time dependence in dV/dx drops out, making dV/dx a function 
of x only (this requires a further redefinition of the time coordinate in the HSS model). 
Generically, dV/dx will depend on both x and Rqr- 

To show that the singularity is a generic feature of these models, we write dV/dx as 

BV 1 Bx 

which, by using the definition of x, can subsequently be written as 

l^,U-R.^, (2.1.) 

oR 3 Xo 

By integrating this expression, we obtain the potential as a function of R and -RgR; 

V{R, Rgk) = - Rgr) -4^ + KRgr), (2.20) 

where A(i?GR) is an arbitrary function of -Rgr- Prom (2.20) and the expressions BV/Bx = 
{R — -Rgr)/3 and 1 + x = x'iR)/Xoi if is clear that for any model in which x'(i?)i? — t- 
(more rigorously, Rx"{R) is integrable) and x(i?) — as i? — ^ oo, the potential will possess 
the singular point BV/Bx oo and V — )• A(i?GR) as x — )• — 1. 

We note that the singularity occurs for models in which F"{R) > at large curvatures, 
and hence is unrelated to the Dolgov-Kawasaki instability [13]. The Dolgov-Kawasaki 
instability corresponds to exponentially growing scalaron modes that are generically present 
in models that satisfy F"{R) < in some dynamically accessible regime. We only consider 
models for which F"{R) > for i? > R^g^. 

2.3 Determination of the Hubble parameter 

In obtaining x and hence R from (2.8) and (2.9), it was assumed that the Hubble parameter 
could be written as H = Hqr + SH, where 5H <^ i^OR is small throughout the evolution 

and hence can be neglected. To check that this assumption is valid, the (0, 0) component of 
the Einstein equations should also be solved for H to check that it does not diverge when 
i? — > oo. 

The (0, 0) component is given by 

18HF"{R) (h + 4Hh) + -3(h + hA F'{R) = (2.21) 
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which, for the AB and HSS models, can be expanded as 



ISH^'iR) [h + AHh) +SH' + ^-s(^H + H^) ^{R) - = (2.22) 

for R > i?vac- Eqn (2.21) is the first integral of the trace of the gravitational field equations, 
and is a second order differential equation for the Hubble parameter. However, due to the 
oscillatory behaviour of i?, it is difficult to solve (2.22) over a significant dynamical range, 
since the divergence in R is extremely sensitive to the initial conditions of H and H. This 
is a manifestation of the singular behaviour of the Ricci scalar. However, we will be able 
to solve (2.21) for the regularized models presented in section 3. 

2.4 Structure of a singularity with F"{R) = for a finite R 

Thus, we have shown that the "Big Boost" weak curvature singularity (2.13,2.14) arises 
whenever RF"{R) is integrable for i? ^ oo with F"{R) > 0. Now, for completeness, let 
us present the structure of an even weaker singularity arising when F"{R) becomes zero at 
some finite value of i?, R = Rg, so that the second stability condition (1.6) is marginally 
violated. For a generic case with F"'{Rs) ^ 0, as follows from Eq. (1.3), this occurs at a 
finite moment of time t = tg when a, H and R remain finite, but R diverges (x \t — t^l"^/^. 
Thus, the scale factor has the following behaviour for t — > tg'- 

a{t) = ao + ai{t - ts) + a2{t - tgf + a2,\t - tsf/"^ + ... . (2.23) 

The metric is C^, but not C^, continuous across this singularity, and there is no unam- 
biguous relation between the coefficients 03 for t < ts and t > ts- 

Therefore, all that was said in section 2.1 regarding the "Big Boost" singularity also 
applies to this weak singularity. Namely, its appearance in solutions of interest for cosmol- 
ogy results in the loss of predictability of future Cauchy evolution. Therefore, if we choose 
to remain inside the scope of F{R) gravity, models where the second of the conditions (1.6) 
is even marginally violated during evolution should be avoided. 

3. Avoiding the weak singularity and solving the problems of F{R) DE 
models 

In the previous section, we have reviewed the existence of the weak "Big Boost" singularity 
which occurs in general at a finite rcdshift in the HSS and AB models, as well as in any 
other F{R) model of present DE for which F"{R) > and RF"{R) is integrable for 
R — )■ 00. Since the problem arises at large curvatures, it is clear that F{R) ~ i? — const 
as — )• 00 is an inappropriate law for viable F{R) DE models.^ Under the less restrictive 
conditions F"{R) > and RF"{R) as R ^ 00, another difficulty, the first point 

*This refers also to combined models of primordial and present DE considered recently in [53, 54] which 
have the same behaviour of F{R) for R —¥ 00. Generically these models possess weak singularities either of 
the type (2.13,2.14), or of the type (2.23), since the second of the stability conditions (1.6) is violated for 
them, too. 
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mentioned in the Introduction, arises; an unlimited growth of the scalaron rest-mass Mg = 
(3F"(i?))~^/^ in the quasi-GR regime (1.7). If Ms exceeds -Mpi, while \R\ remains less 
than -M|j, particles (scalarons) should collapse to black holes from a naive physical point 
of view. More formally, this means that loop quantum-gravitational corrections to the tree 
action (1.1) become dominant, so it may not be used further in an effective quasi-classical 
theory. Thus, such models of present DE are in general not compatible with the standard 
early Universe cosmology including the correct BBN and recombination. 

Since the origin of these two difficulties is that F"{oo) = 0, they can both be cured by a 
very simple change in the HSS and AB models. Let us add an additional term, quadratic in 
the Ricci scalar, to their F{R) functions, so that F"{oo) becomes non-zero. In particular, 
we consider the following functions, 

Fhss = Fhss{R) + -^, Fj,b = Fab{R) + -^, (3.1) 

where M is a mass scale coinciding with the scalaron rest-mass whenever low curvature 
modifications to GR can be neglected. There is no "Big Boost" singularity in solutions of 
the i?^-corrected HSS and AB models (3.1). 

If we now return to the linearized analysis for these new models, in which R = Rqr + 
SR, then it follows from Eq. (2.5) that the quadratic term has two important effects on the 
dynamics of SRosc - the oscillating component of 6R. The first is that it introduces an upper 
bound on the mass of the scalaron, Mg < M, and hence limits the frequency of oscillations 
of R (as was noted in [44]). The second effect is to moderate the amplitude growth of SRqsc, 
which now goes oc a~^/^. Specifically, going to the past, 6Rosc/Rgr decreases throughout 
the matter and radiation eras beginning from the moment when the R^ correction in 
Eq. (3.1) becomes larger than the non-GR term decaying with the growth of R. Hence 
SRosc remains a small perturbation throughout the cosmological evolution, subject to it 
being small at present. Vice versa, going to the future time direction, since SRosc/ Rgr 
grows until recently for the models (3.1), it remains an open problem to explain its very 
small initial amplitude in the early Universe. It is clear that the scalaron overabundance 
problem noted in the Introduction has not yet been solved. 

Let us now discuss possible values of the parameter M in Eq. (3.1). It should be 
sufficiently large in order to satisfy the viability conditions presented in the Introduction. 
In particular, it may not be ~ y/Rvac ~ 10~^^ eV as was considered in [41, 55], and even 
the values discussed in [56, 57] are not high enough to solve the overabundance problem. 
A non-oscillating part of SR induced by the R^ correction to F{R) becomes important 
when H(t) ~ M. As a consequence of this, the lower limit M > 10^^'^ eV which follows 
from the most recent laboratory Cavendish-type experiment [58] is already sufficient for 
this correction to GR to be negligible both during BBN and in the center of neutron stars.^ 

^Still, as we shall see from section 4, non-GR terms in Eq. (3.1) might become important should the 
trace T of the matter energy-momentum tensor change sign and become positive inside neutron stars, as 
happens in idealized p — const solutions considered in [59, 57). However, it is argued in [60] (see also [61]) 
that this does not occur inside realistic neutron stars. 
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However, taking M close to this lower limit is incompatible with the existence of any 
kind of inflation (not specifically driven by F(R) gravity) in the early Universe, which is 
required to solve many other cosmological problems. Indeed, if M <C i^inf where i^inf is the 
Hubble parameter during the last part of inflation, quantum fluctuations of the scalaron, 
including long- wave ones, are generated during inflation that generically results in a large 
value of SRosc after the end of inflation, and can even lead to the existence of a second 
stage of inflation driven by the scalaron itself [62]. Of course, the scalaron is not stable and 
decays into pairs of particles and antiparticles of all non-conformal quantum fields (this 
is a particular case of the effect of particle creation in gravitational fields). However, this 
process is sufficiently slow. Even in the pure R + i?^/6M^ model which does not describe 
the present DE, the characteristic scalaron decay time is r ~ AfpjM"^ [8, 27, 28, 31]. 
Thus, one would need M S> 10^ GcV to ensure the scalaron decay by the moment when 
~ 1 s~^ to avoid any problems with BBN. However, we shall see in the next section that 
a non-trivial low-R structure of F{R) models describing present DE results in slowdown 
of the scalaron decay after the end of inflation driven by the scalaron itself. Since the 
aim of the present paper is to find at least one F{R) model of present DE satisfying all 
five viability conditions formulated in the Introduction, we choose the value of N which is 
sufficient to solve the scalaron overabundance problem. Namely, we take M either larger 
than i^iuf in case inflation is produced by some other scalar field, or ^ 3.7 x 10^''^ {50 /N) 
GeV if the scalaron plays the role of an infiaton, as was noted in the Introduction. In the 
former case, scalarons are practically not generated during inflation, so 6Rosc is zero after 
its end. In the latter case, reheating after inflation appears to be very non-trivial; it is 
studied in section 4. 

Let us now return to the evolution of the i?^-corrcctcd HSS and AB models at recent 
redshifts. Although the above reasoning suggests that the perturbative analysis is sufficient, 
it is also interesting to investigate the effect of the quadratic term on the scalaron potential. 
Performing the same steps as in section 2, the potentials are exhibited in fig.2. We see 
that the singular point at which V{x) — >■ const, dV/dx — >■ — oo is no longer present, and we 
observe a regular potential for all values of x of physical interest. 

In fig. 3, we have confirmed numerically that the introduction of the R? term bounds 
the oscillations, which now grow to the past at a slower rate than .Rgr in the AB model. 
To obtain these curves we have numerically evolved the trace of the gravitational field 
equations, using the i?^-corrected AB model and taking 5 = e/M^ = 4 x 10~^. We 
have evolved equation (1.3) through the matter era, using a dimensionless time coordinate 

1/2 . ... 

t = Rykct and using random initial conditions for R and dR/dt. In the AB and HSS models 
without the R^ term, it was found that R could only be evolved over very short timescales 
(evolving backwards over the matter era), and only if we fine tuned initial conditions to 
approximately R ~ -Rgr- Now, we can choose random initial conditions for Ri, Ri and Hi, 
Hi and evolve backwards. In all cases we observe that the oscillations of R decay to the 
past, and no singularity is present. Unfortunately, wc cannot directly compare the results 
obtained here with the corresponding functions 5R and 5H in the uncorrected AB and 
HSS models. This is due to the fact that the uncorrected models will generically evolve to 
a singularity over much shorter timescales than presented here (see ref.[46] for a discussion 
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Figure 2: The potential for the i?^-corrected AB model (top) and HSS model (bottom). We have 
taken & = 4, i?GR = 10 x i?vac, 5 = e/M^ = 1.5 x 10"^ and ai = for the AB model, and n = 2, 
c = 1, i?GR = 10 X i?vac, 5 = 1.5 X 10~^ and a2 — for the HSS model. We see that V{x) — oo as 
X — >■ —CO, and hence the singular behaviour has been removed from both models. Specifically there 
is no singularity at x —> —1, since both V{x) and dV/dx are now regular at this point. 

of the difficulties associated with numerically modelling the original AB and HSS models). 

Summarizing, we have found a way to cure all three problems of the HSS and AB 
models of present, low-curvature, DE which does not destroy correct inflation, BBN and 
other advantages of the early Universe cosmology. The approach consists of changing the 
behaviour of F{R) at R ^ i?vac according to Eq. (3.1) with a very large value of the free 
parameter M which should either exceed the scale of inflation, or be equal to the concrete 
value needed for scalaron driven inflation in F{R) gravity. Still this not the end of the 
story, one loophole remains which requires further correction of the functions (3.1), now in 
the range R < i?vac including negative values of R. This final step in constructing a viable 
F(R) DE model will be made in section 4. 

Let us finish this section with a comment on possible alternatives to the large- be- 
haviour (3.1). If instead a more general term M^^^^i?™, m > 1, is added to Fhss or Fab, 
it can be checked that Mg, while growing with R, never exceeds it, so it may not become 
larger than Mpi if R is less than Mpj. Thus, there is no problem with the unlimited growth 
of Mg. However, instead we face a new difficulty: a possibility of the formation of a new 
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R_{vac)'^(l/2)t 



R_{vac)'^(l/2)t 



Figure 3: The fractional difference 5i?/i?GR = {R— i?GR)/^GR for the i?^-corrected AB model, 
obtained by numerically evolving the trace of the gravitational field equations using S = e/M^ = 
4 X lO^'*, and taking i?GR — 4/3t^ (that is, we evolve the Ricci scalar assuming matter domination). 
We take random initial conditions for R and dR/dt, and use a dimensionless time coordinate 

1/2 

t = Rvact, where t is the cosmological time, (a) and (b) differ only by the initial conditions placed 
on R and dR/dt; in (a) we have perturbed R significantly away from its General Relativistic value. 
We observe that the oscillatory component of the Ricci scalar decays to the past relative to -Rqr, 
as predicted in the text. 



space-like curvature singularity with 

a{t) oc {ts -ty , q 



{m - l){2m - 1) 



m 



(3.2) 



during evolution to the future which destroys all subsequent evolution [63, 64, 65]. For 
m > 2, this singularity is of the "Big Rip" type, and the scale factor a{t) becomes infinite 
at a finite moment of time. For 1 < m < 2, the scale factor becomes zero at t = ts- In both 
cases, it remains an open problem if the singularity (3.2) can be avoided in generic future 
Cauchy evolution after inflation. So, at present the large- i2 behaviour (3.1) of F{R), up to 
logarithmic in R corrections, seems to be the only one free from dangerous pathologies. 



4. Inflation and late time acceleration from one F{R) function 

As has been shown in the previous sections, to construct a viable model of present DE in 
the scope of F{R) gravity which satisfies all the viability conditions and does not destroy 
previous successes of the early Universe cosmology, one has not only to choose the correct 
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non-GR structure of F{R) at low curvatures R ~ Rvac, like that in the HSS or AB models, 
but also to modify high-i? behaviour of F(R). Moreover, it appears that the only high-i? 
behaviour which does not lead to problems with new singularities is just that was originally 
proposed for the scalaron driven R^ inflation in the early Universe. Therefore, it is natural 
to consider combined F{R) models which describe both primordial and present DE using 
one F{R) function, albeit one containing two greatly different characteristic mass scales. 

It presents no problem to find a function F{R) for which the equation (1.4) defining 
de Sitter solutions has two or more roots. The real issues are, first, to ensure that the 
inflationary de Sitter solution is metastable, slow-roll and leads to the correct spectrum of 
primordial perturbations and, second, that there exist a sufficiently effective mechanism of 
reheating after inflation which transfers energy from scalarons into ordinary matter and 
radiation and heats them to a high temperature long before the beginning of a second 
de Sitter stage which wc observe now as the present acceleration of the Universe. This 
transition between the two accelerating epochs in F{R) gravity must be carefully analyzed 
and the absence of singular points or instabilities in the cosmological evolution has to be 
proved. 

4.1 New problem 

When wc begin to consider post-inflationary evolution in F(R) gravity a new problem 
immediately arises, not only in the combined case when inflation is scalaron-driven but 
also when inflation is produced by a minimally coupled scalar field ip, which requires 
further generalization of even the corrected models (3.1). Namely, while R is positive 
during both inflation and the recent evolution of the Universe, it becomes negative (and 
large) during each post-inflationary oscillation of ip. In particular, T = 3P — p = -i/'^ > at 
the moment when V{ip) = 0. It will be shown below that the same occurs after scalaron- 
driven inflation. Thus, the range of R in the models (3.1), previously used for R > Ro only, 
has to be extended to negative values, and the stability conditions (1.6) should be satisfied 
for negative R, too, at least up to values R ~ — which occur during post-inflationary 
evolution. 

However, this is not possible to achieve without a further change of these models. 
In the case of the HSS model, both original and the ii^-corrected one, flrst, one has to 
assume additionally that n is an integer to avoid non-analytical behaviour at i? = (note 
that there is no such problem in the variant of this model introduced in [44]). A more 
serious problem is the appearance of the weak singularity (2.23) at the points R/Rya^ = 
± ((2n - l)/c(2n + l))^/2n ^j^g^^g F"{R) = 0. As shown in section 2.4, it is not possible to 
predict deterministic Cauchy evolution through it in the generic case (when 03 7^ 0). 

In the case of the AB model, its function Fab{R) is analytic and satisfles the conditions 
(1.6) for all R. However, F^^(— 00) = and it approaches this limit exponentially fast for 
low values of \R\: F'{R) ~ exp[2(i?/e — b)] for < 0, \R\ S> e. As a consequence of this, 
the i?^-corrected AB model (3.1) acquires a point at which F'{Rq) = for Rq ~ —Rvac- As 
explained in the Introduction, such points should be avoided because of the formation of a 
generic anisotropic curvature singularity [10, 11]. Since this occurs at very low curvatures 
in this case, this point is dynamically accessible during the post-inflationary phase of the 
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Universe evolution, where R oscillates around the vacuum state. Hence, the ii^-corrected 
AB model is not viable either. 

Moreover, integrating the viability condition F"{R) > over the interval {—R, R) with 
R <^ M^, we obtain that any viable F(R) DE model should have a non-zero ^-factor 

F'{R) - F'{-R) 



g 



(4.1) 



2F'{R) 

Physically this means that the value of the effective background Newton gravitational 
constant G^s = G/F'{R) in the quasi-GR regime (1.7) is larger for < than for i? > 0. 
It follows from the stability conditions (1.6) that the (7-factor always lies in the range 
< 5 < 1/2. 



0.40 



0.35 



0.30 



10 



20 



30 



40 



50 



Figure 4: The allowed parameter range of g and b. For any choice of h and g above the curve, the 
model (4.2) has a stable de Sitter vacuum state. The curve is obtained by calculating the values of 
g and h for which the functions Q{R) = RF'{R) — 2F{R) and dQ/dR are both zero. 



4.2 Resolution of the problem and the improved AB model 

This new problem can be solved and the problematic point F'{R) = for small can be 
avoided if we consider the improved, ^-extended i?^-corrected AB model, more concisely, 
the gR'^-AB model. 



F{R) = {l-g)R + gelog 



cosh(i?/e-6) 



cosh b 

where we have introduced the new dimensionless parameter - 17- factor, < 5 < 1/2. Note 
that g = 1/2 corresponds to the i?^-corrected AB model, and (7 = to the R? inflationary 
model with F{R) = R + R'^/6M'^ [44] which does not present DE. The function (4.2) 
corresponds to an interpolation between two different gravitational constants, as F{R) ~ R 
for i?vac < R< and F{R) ~ (1 - 2g)R for -M^ < R < e, with a step at ~ be. 

4.3 de Sitter attractors 

Like the AB and i?^-corrected AB models, this new model can be expanded as F{R) ~ 
R - i?vac/2 + xiR) for i?vac < -R < M^, where now i?vac = 2(7(6 + log[2 cosh 5])eAB and 
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xiR) = 5'ee^''e~^^/^ + i?^/6M^. It has stable Minkowski and de-Sitter vacuum states for an 
appropriate choice of b and g; in fig.4 we have exhibited the aUowed parameter range. To 
obtain this curve we have used the fact that for particular values of b and g, the function 
Q{R) = RF'{R) — 2F{R) will possess three zeros corresponding to vacuum states of the 
model (Minkowski space, and two de Sitter vacua). The values of g and b below the curve 
yield only one zero in Q{R) (Minkowski space), and those above will yield three. The 
curve corresponds to the limiting case where Q{R) possesses a double zero, that is when 
Q{R) = and Q'{R) = 0. Hence we find that for any parameter choices of 5 and g above 
the curve, there exists a stable de Sitter vacuum state. As we increase 6, the allowed range 
of g increases, and as 6 ^ oo, g must satisfy g > 0.25. We also note that as 5 ^ 0.5, the 
condition b > 1.6 is required for a de Sitter state. Finally, the model satisfies F'{R) > for 
R > -3M2(1 - 2g) and F"{R) > for all R, and hence possesses no known instabilities 
for R > -3M2(1 - 2g). 

4.4 Slow-roll inflation 

We now study the evolution of the model (4.2), starting from a high curvature slow roll 
epoch and evolving forwards in time. We assume that the classical evolution of this model 
begins with Hi ^ Mpi and \H\i ~ ^ Hf. We analyze the (0,0) component of the 
Einstein field equations 

18HF"{R)H + 72F"{R)H^H + - 3 (^H + H^^ F'{R) = 0, (4.3) 

in order to study the behaviour of the Hubble parameter. In the regime < R < Mpj, 
the function (4.2) has the form F{R) « R+R^ /6M'^ - R^i,c/2 + geex.p[-2{R/e-b)], and the 
last two terms are negligible (in particular, the low curvature correction is exponentially 
suppressed). Therefore it is an excellent approximation to use F{R) = R + R'^ /QA'P in 
(4.3), which gives 

H--^—L + SHH+—H = 0. (4.4) 

The slow-roll evolution of (4.4) (and generalizations thereof) has been studied extensively 
[33, 28, 29, 30], and it has been shown that for > \H\, H\H\ > \H\, the Hubble 
parameter has the following form 

H{t) ^Hi--^, (4.5) 

where we have taken t = at the beginning of the evolution and Hi ^ Mpi. Slow roll ends 
when H ~ M. Following this, we have an epoch in which the both the Hubble parameter 
and Ricci scalar oscillate with high frequency around a stable ground state. In the following 
section, we consider the oscillations of the Ricci scalar in detail, and how they may reheat 
the Universe. 

Let us show now that all other scalaron-driven infiationary models in F{R) gravity 
which produce primordial spectra with other values of and r are, in a sense, close to the 
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model. Indeed, there may be two kinds of inflationary models. In the first one which 
is the analogue of chaotic inflation [66] in GR, inflation occurs over a wide range of R. 
Then it follows from Eq. (4.3) that we need F{R) ^ R^A{R) for i? ^ oo with A{R) being 
a slowly varying function of R, namely 



\A'{R)\ < 



A{R) 



\A"{R)\ < 



A{R) 



(4.6) 



R ' R'^ ' 

for this kind of inflation to take place. Thus, these models are indeed close to the R? one. 

In the second case, inflation occurs around some fixed root i? = i?i of Eq. (1.4) - an 
analogue of the "new" inflationary model [67, 68] in GR. Then, taking into account that 
the inequality (1.5) should be satisfied only marginally for metastability of the de Sitter 
solution, we get 



Ri 



F"{Ri) 



Rl 



(4.7) 



Thus, these models are close to the R^ one near the point R = Ri. 

Finally in this section, we exhibit the potential of the scalaron in the Einstein frame 
in fig.5. We observe the metastable de-Sitter point, and the stable Minkowski vacuum. 

4.5 Reheating 

We now study the reheating epoch for this class of models (a discussion of reheating 
mechanisms for modified gravity models can be found in, for example, [27, 28, 31]). For 
the R^ inflationary model, immediately following slow roll we have ~ |ij| ~ M^, and 
the Ricci scalar undergoes damped harmonic oscillations around -R = with frequency 
uj = M and \R\ <^ M^. In this section, we consider the evolution of H for the model in 
question, taking initially ~ M^. We find behaviour that is dramatically different 

to that of the R^ model, suggesting that the modifications to General Relativity at ~ 0, 
which were introduced to induce late time acceleration, also have a significant effect on 
the dynamics of the Ricci scalar in this early epoch of the Universe. In any unified model 
of inflation and present DE there must be an efficient reheating mechanism, otherwise the 
Universe will relax to its late time de Sitter attractor without first undergoing epochs of 
radiation and matter domination. 

To study reheating for the model (4.2) we again consider the (0,0) field equation, which 
can be written in terms of dimensionless parameters H = H/M, R = R/M'^, S = e/M^ <^ 1 
and dimensionless time t = Mt, 



HH" - + 3^2^' + ^2 -1{H' + H^) tanh 



12 



log 



cosh{R/S - b) 
cosh(6) 



+ ^sech^ 





R 
J 



{HH" + AH'^H') = 0, 



(4.; 



where primes now denote derivatives with respect to i. Initially, we will solve equation (4.8), 
neglecting any effects due to gravitational particle production (that is, we will neglect the 
backreaction of created particles on the dynamics of H{t)). The issue of backreaction will 
be tackled in section 4.5.2. 
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Figure 5: The Einstein frame potential for tire gR^-KQ model (solid) and the original AB model 
(dashed). To obtain the potential we have used g — 0.47, 5 ~ e/AP = 10^"* and h — 2. We observe 
two vacuum states; one late time de Sitter attractor and a flat Minkowski vacuum state. Slow-roll 
occurs for (f) ^ Mpi. We note that in the original model, the scalar field satisfies < and the 
singularity corresponds to the point (/)—>■ 0. This singularity is removed in the gR^-KQ model by 
mapping the singular point i? — > oo in the Jordan frame to — > oo in the Einstein frame. 



4.5.1 Evolution of H{t) without backreaction 

To begin, we solve (4.8) numerically. We are considering the epoch immediately following 
slow roll, so we take as initial conditions H[ = —1 and H\ = \ and evolve in t (we will 
begin our evolution at t = 0). The results are exhibited in fig. 7 for the Hubble parameter 
and fig. 8 for the absolute value of the Ricci scalar, where we have taken 5 = 1.5 x 10~^, 
6 = 3 and g = 0.45. We note that for realistic choices of M and e, 5 will be many orders 
of magnitude smaller, however numerically we are restricted to (5 ~ 10^^. As can be seen 
in fig. 7, H evolves through a number of distinct regimes, which we discuss below. 

Initially, ~ 1, so we can use F{R) k, R+E? /QM"^ , and the Hubble parameter behaves 
as in the B?' infiationary model (that is, H rapidly decreases io H ^ 1 over timescales 
t ~ 0(1)). This region corresponds to t < 1 in figs. 7,8. In this regime R ^ S, and hence 
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the value of 6 will have no significant effect on the dynamics of H and R. 

After this (very short) initial period, we observe a significant change in the gradient 
of H. In this regime, we have H'^ :$> S, \H'\ :$> 6 hut R = 6H' + 12H'^ > 6 (that is, 
H' —2H^.) The exact value of R will depend on the specific function F{R) being 
considered, and corresponds to its value when the large curvature contribution to the 
scalaron mass is approximately equal to the low curvature terms; for the model constructed 
here when 1/3M^ ~ {g/e)scch'^[R/€ - b]. 

To analytically model the behaviour of H in this regime, we look for a solution to (4.8) 
of the form H{i) = Ho{i) + dHi{i) + OiS"^). In doing so we obtain 



6 {HoH'^-i^ + 3H^Hl,+ 



2 £y/ , (1 ^^H^-^ {H'o + Hi) tanh \6H[ + 24HoHi - b 



+ S^fsech 



2 ■ - u u . 2 "2 
6H[ + 24HoHi - b] {HqH'^ + AHIh'^ 



+ 35(5sech2 m[ + 2AHqHx - b {HqH'( + H^H'^ + SHqHiH'o + 4HIh[) = 0, 



(4.9) 



where we have anticipated that Rq = 6-ffg + 12f^Q = at zeroth order in the tanh[i?/(5 — b] 
and sech^[i?/(5 — b] terms (we will confirm that this assumption is valid below). At order 
0(1) and 0{d) we arrive at 



0(1) HoH'f; + 4HlH'o = 0, (4.10) 

+ 2^ 



cosh^ 



6H[ + 2AHqHi - b 

0{5) HoH'l + HiH'^ + moHiH'^ + 4HIh[ = ^ ^ x (4.11) 



- i?oi?o - 3i?o'^o - ^^^^0 + ?(^o + ^o) tanh 6H[ + 2iHoHi - b 



and hence at 0{1), Hoit) = l/(2t + ao), where ckq is an integration constant. At zeroth 
order we have a{i) oc i^/^, H{i) ~ l/2i and R = 0. This is in agreement with our numerical 
results; in fig. 6 we have exhibited H as calculated numerically and H = l/2t logarithmically 
as a function of log[f], over the range log[t] = (log[l], log[180]), and we see a close agreement 
between our analytic and numerical results. The number of scale factor e-folds during this 
period is given by A''! '-^ — ln(l — 2g). 

Reheating in this class of models will begin at this stage in the evolution of H, via 
gravitational particle production due to the abrupt change in the Ricci scalar (in our 
numerical calculation this occurs at t ~ 2, as shown in fig.8). Prior to this, the Ricci 
scalar is approximately in a de Sitter phase, which then abruptly evolves to a radiation-like 
Universe with i? ~ 0. The main contribution to the number of particles and antiparticles of 
all non-conformal matter quantum fields created after the end of inflation occurs during this 
(almost) discontinuous change in R. Using either the expression for the rate of creation of 
the density of scalar particles and antiparticles n in the massless limit in a FRW background 
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Figure 6: exhibits H (solid), obtained by numerically solving (4.8), and \/2t (dashed) plotted 
logarithmically against t. We see that over the range t — (5, 180), H ^ l/2t, in agreement with our 
analytic result. 
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Figure 7: The behaviour of H using the parameters 5 — 1.5 x 10^^, h = g ~ 0.45, and neglecting 
backreaction. We observe a number of distinct regimes; for < t < 2 the Hubble parameter is 
exiting the slow roll regime. Following this is an epoch for which H ^ l/2t. The abrupt change in 
Hati^ 180 corresponds to a spike in the evolution of R. Following this H ^ S^^^ and a{t) ^ const. 
Once H spikes once more a,t t ^ 5000, the Hubble parameter has completed one full oscillation. 



[69, 28] 

which for ^ = is vahd for gravitons (although the rate should be doubled in this case 
due to two polarization states) and the longitudinal component of vector bosons, too,^ or 



^Strictly speaking, Eq. (4.12) is derived in the limit < i^^^i?'"'. 
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Figure 8: The behaviour of the Ricci scalar R using the parameters 5 = 1.5 x 10~^, b = 3, g = 0.45, 
and neglecting backreaction. We observe the Ricci scalar evolving from a slow-roll de Sitter state 
to R> 6, where the low curvature corrections to GR dominate. R is typically ~ d throughout this 
epoch, although it periodically exhibits spikes. 

results of the papers [52, 70], we arrive to the following estimate of the total number and 
energy density of created ultrarelativistic particles: 

n ~ xH^ , prad = xHf , (4.13) 

where H^- is the value of the Hubble parameter at which the discontinuous change in R 
occurs (for our model H,. ~ M), x is a dimensionless parameter of order x ~ 10~^ and 
is the value of the scale factor at this point. The energy density of created particles will 
be initially subdominant, but may grow to have a significant backreaction effect on the 
dynamics. For the remainder of this section we neglect the effect of prad on the evolution 
of H and consider backreaction in the following section. 

Over the regime H = l/2t + 0{5) discussed above, the Ricci scalar is given by 
i? ~ 6(5 (^H[ + AHqHi^ , which decreases from i? > 5 to i? ~ 0, at which time the 

(1/5) sech^(i?/J— 6) term in (4.8) no longer dominates. At this point, we can use sech^(i?/5— 
6) ~ and tanh{R/6 — 5) ~ —1, in which case (4.8) becomes 

HH" - + 3H^H' + ~ - 0, (4.14) 

which has solution H ~ Aosm'^{^/T^^2^i/2) (that is, H oscillates on timescales of order 
t ~ 1, as in the inflationary model). This corresponds to the spike in 7? at t ~ 200, with 
amplitude R > —1. However, unlike in the R^ inflationary model, R only completes one 
half of an oscillation before we once again enter a regime where the (1/6) sec]i^{R/6 - b) 
term dominates. 

Following the spike in R, H does not decay like H cx t~^. Instead, H^,H' < 6, and the 
scale factor a(t) is approximately constant over this regime. This is confirmed numerically 
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in fig. 9, where we observe a plateau in the evolution of the scale factor. During this period, 
the number e-folds will be N2 ~ 0, and H and R will grow until R > S. 

Following this, R will once again produce a spike (with amplitude R < 1), completing 
one full oscillation of H. Numerically we can only observe one complete oscillation before 
R reaches its final state R > S, however this is simply a consequence of being unable to 
choose a sufficiently small value of S. Finally, we note that if we average H over an even 
number of oscillations, we find (H) = l/3t, as one might expect for a period of kination 
(domination of a masslcss scalar field). 

We can summarize the dynamics of H as follows; the Hubble parameter periodically 
undergoes (almost) discontinuous jumps between periods when H = l/2i and H ~ (5^/^. 
The duration of these periods in terms of Int is determined by the non-zero ^-factor (4.1) 
and is equal to ~ —2 ln(l — 2g) and — ln(l — 2g) respectively. It is clear that such behaviour 
of H and R is markedly different to the standard reheating dynamics, and we expect 
that the low curvature modifications to General Relativity will leave unique observational 
imprints. In particular, our reheating mechanism is less efficient than the pure F{R) = 
R + R^/6M'^ model, and we expect a significantly lower reheat temperature. Additionally, 
the average expansion rate after inflation is slower than a{t) oc t^/^, and therefore the 
number of e-foldings, A'^, should be ~ 70, larger than in standard inflationary models and 
the same that would occur if all inflation proceeded at H = Mpi . Indeed, by comparing 
the energy density of created particles (4.13) to H^, we see that they become equal at 

t = treh ~ x-^/^M-^Mli ~ 10-^^ s (4.15) 

after the end of inflation (assuming the value M ~ 3 x 10^'^ GeV needed to fit the amplitude 
of observed curvature fluctuations with N = 70). If g = 0.45, then about 7 complete 
non-linear oscillations of R have occurred by this moment. If the radiation component 
has been already thermalized at t = treh due to interactions between particles, then its 
temperature is r(treh) ~ 10^ GeV which is sufficiently large (however, in principle, it may 
thermalize significantly earlier while being a sub-dominant component). At t = tj-ch, the 
comoving scale which was equal to the Hubble radius at the end of inflation is given by 
M-^Mtreh^/^ = x-V2m-2Mpi = x^/^Mp^^^tl^l. Therefore, up to a factor of a few, it 
coincides with the comoving scale which is equal to the Planck length at the Planck time 
in a universe which is radiation-dominated at subsequent times. As a result, irrespective 
of the fact that if the thermodynamic equilibrium in the radiation component is reached 
before or after tj-chj this scale coincides with the characteristic thermal length of present 
CMB photons with temperature = 2.725 K, up to a purely numerical factor depending 
mainly on an effective number of species at the moment when the equilibrum has been 
achieved. This explains why N ^ 70 for our model. In the next subsection, we consider 
backreaction of created particles of a FRW background numerically and in more detail. 

Hence, due to this change in N, the index (slope) of the power spectrum of primordial 
scalar (density) perturbations Hg is slightly higher, Us = 1 — 2/N ~ 0.97, in our model which 
combines both inflation and present DE using one F{R) function (4.2), as compared to 
the R + R'^/QM'^ model describing inflation only. However, this distinctive and observable 
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prediction is degenerate, since it may be changed by introducing loop corrections to the 
large curvature E? term of another purely inflationary model, making it of the type (4.6). 
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Figure 9: The behaviour of a(i), normalized such that a{ii) = 2 x 10 ^ . We observe a plateau in 
the evolution of the scale factor, as predicted in the text. 



4.5.2 Effect of backreaction 

As pointed out in the previous section, at the beginning of the reheating epoch R undergoes 
an almost discontinuous change from a de Sitter phase to R> 5, and particle production 
occurs mainly during this period. To take into account the backreaction of these particles 
on the evolution of the Hubble parameter, we must now solve the equation 



{H 



i\2 



(1-5) 



HH" - + 3H^H' + ^^^H^ -^(H' + H^) tanh 



-\ lOE 

12 ^ 



cos\i{R/5 - b) 
cosh(6) 



H — ^sech^ 





R 



(HH" + 4i^2i^'^ 



Prad 



(4.16) 

xM^at 



For future convenience we define the dimensionless parameter /5rad = Pra.d/M'^. Initially, the 
energy density prad of produced particles will be subdominant to the background evolution 
of H{t), driven by the scalaron. Numerically, we will take the radiation component to be 
Pi/M^M^i ~ IQ-'^H?, and xM'^/6M^^ = 2 x 10"'* (one must be careful to preserve the 
hierarchy between H, prad and 5, as when prad is of order 5, our solution will cease to be 
physically relevant.) 

We begin our numerical evolution at the point where the radiation energy density is 
produced, so M^pi/6M^^ = xM^/6M^^ = 2 x 10"^, and use initial conditions 2H^ = -H[. 
Taking ,5 = 1.5 x 10"^, 6 = 3 and 5 = 0.45, we present H"^ /2 and M'^p/6M^^ in fig.lO. 
We note that the observed behaviour of H is significantly different to that of the previous 
section, however our conclusions will remain essentially the same. As before, we observe 
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Figure 11: The Hubble parameter during one complete oscillation (solid). We have also exhibited 
the two approximate analytic solutions as calculated in the text; H = l/2t (dashed) and H = 
l/{2i+a) (dotted). 



one complete oscillation of H before p ^ 5. We now discuss the various regimes over the 
course of this oscillation. 

Initially, H behaves in a very similar manner to the previous section; the radiation is 
subdominant and H = l/2t, p cx t~'^. During this time R evolves from i? > 5 to ~ 0. 
Then, at the point i? ~ 0, the Ricci scalar undergoes half an oscillation, spiking at i? > — 1 
(again, as before). However, H dies not oscillate to -ff = as in the previous section, but 
rather H'^ oscillates around M'^ p/3Mp^. We see that by incorporating p into the dynamics, 
H now satisfies H ^ throughout the reheating epoch. 

All that remains is to consider the behaviour of H and a{t) over the second half of the 
oscillation of H. Following the (almost) discontinuous change in at i ~ 10, the Hubble 
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parameter satisfies H ^ \f5 and H' ^ 5, however R> 5. By using the same perturbative 
analysis as before, we must conclude that H is given by H = l/(2t + a), where a is an 
integration constant. The value of this constant can be obtained numerically, and depends 
on the ratio of Hf and M'^p/Mp,. For realistic values of M and Mpi we expect that a will 
be much greater than t, and hence H ~ and a{t) Ri const during this period. The 
number of e-foldings over this regime will be N2 ~ 0, as before. During this time R grows 
from i? ~ to i? > (5, at which point we observe another spike in R, thus completing 
one full oscillation of H. In fig. 11 we have exhibited the Hubble parameter and the two 
analytic approximations l/2t and l/{2t + a) derived above; we see that our results are in 
agreement with the numerical solution. 

Although the evolution of H is significantly modified when we incorporate backreac- 
tion, our conclusions remain essentially unchanged. H now undergoes periodic oscillations 
around ~ M'^p/Mp,; the first half of this oscillation is characterized by the behaviour 
H = l/2t, whereas its behaviour over the second half is given hy H = l/{2t+a) « 1/a <C 1. 
As in the previous section, we find that the reheating mechanism is less efficient than the 
standard R + R^/6M^ model, and the averaged Hubble parameter will evolve at a slower 
rate than l/2t, so our conclusions regarding the total number of e-foldings and the conse- 
quent increase in the value of ng will persist. 

4.6 Cosmological evolution 

Following the reheating epoch of the Universe, we expect that the model (4.2) reproduces 
the standard cosmology, that is it evolves from a radiation dominated epoch to matter 
domination, with the final state of the Universe being the de Sitter vacuum. To see that 
the standard cosmology is reproduced, we rearrange the (i,j) and (0,0) gravitational field 
equations, assuming that the energy-momentum tensor is comprised of matter and radia- 
tion fluids only, obtaining 



3^2 ^ Pm + P r + Pf{H,H,H,R) 



Ml, 



2H + 3H 



2 Pra + Pr + PF{H,H,H,R) 



Ml, 



(4.17) 
(4.18) 



where the subscripts m and r represent the matter and radiation components, and F 
denotes effective density and pressure terms due to the F{R) modified gravity function (we 
have made it clear that pp and Pp depend on H and its derivatives, to stress that we are 
dealing with a system of fourth order differential equations), pp and Pp are given by 



PF 

Ml, 



3H 



3M2 



+ ^sech2(i?/e - b) 



R 



R' 



12M2 



+ (ij + i?^)A + |iog(coshfe) 



- I [e log cosh(i?/e - b) - R] + 'ig{H + H^) (tanh(i?/e - 6) - 1) , (4.19) 
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^''^%anh(i?/e-6)sech2(i?/e-6) (4.20) 



Ml, e2 



+ ( + -^sec\i^{R/e -b)]{R + 2HR) - g{H + ZH^){iwh{R/e - 6) - 1) 



R — e log 



cosh {R/e — b) 
cosh 6 



R R 



2 



3M2 12M2 ■ 



It is clear that pp and Pp will act as dark energy components, and by solving the gravita- 
tional field equations we can obtain the equation of state parameter w = Pp/pF- 

To begin, we note that previous studies [44] have shown that the Ricci scalar can 
be written as R = i?GR + ^Rosc + ^Rind during the cosmological evolution of models 
such as (4.2), where SRosc is the oscillatory, scalaron component, and SRind is given by 
6Rind - DF'{Rgr) + RgrF'{Rgr) - 2F{Rgr) + Rgr- The scalaron oscillations SR^sc 
have been discussed in previous sections; they are well behaved and regular for the model 
(4.2), and we can assume that the energy density of these oscillations has decayed during 
reheating and can be neglected. SRind also satisfies SRind ^ ^GR for R > i?vac, and hence 
as an excellent approximation we may simply use R ~ Rgr and H ~ Hgr in p-p and pp. 

Typically during the matter and radiation era's, we have e ^ -Rgr ^ M"^, in which 
case (4.19) and (4.20) can be written as 



PF 

Ml, 



2 



b + log(e^ + e-^) 



+ 



M2 



{H + H^)R - HR 



12 



(4.21) 



_g-2(i?/6-6) 



Pf 
Ml, 



9± 

' 2 L 



b + log(e^ + e-^ 

^g-2(ii/e-6) 



+ 



12M2 



A{R + 2HR) + R^ - A{H + 3H'^)R (4.22) 



2 



2g{H + 3H^) + ^{R + 2HR) - ^{R? 



Since last two terms in (4.21) and (4.22) are suppressed by factors of l/M^ and er'^W^-^) 
respectively, they are completely subdominant, and throughout the matter and radiation 
era's we have w ~ — 1. This simply reflects the fact that for R S> -Rvao the mass of the 
scalaron Afscai is very large, and as an effective field theory below energy scales E ^ Mgcai 
the model (4.2) reduces to GR with a cosmological constant. 

However, at late times, when Rgr ~ C(-Rvac)) the mass of the scalaron is small, and 
we expect that there may be significant deviations from w = —\. The change in w is due 
to the tanh(i?/e - 6) - 1, sech2(i?/e - b) and log[cosh(i?/e -b)\-R terms in (4.19, 4.20), 
which are no longer exponentially suppressed. This behaviour is exhibited in fig. 12, where 
w is shown as a function of rcdshift. w{z) was calculated by solving the full gravitational 
field equations numerically, starting at z = 4 and evolving to z = 0. As initial conditions. 
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we have assumed that H and its derivatives are exactly their GR counterparts at z = 4, 
and f^m = 0.3, and have taken 6 = 10~^ . With this choice, we observe that w oscillates 
around —1 in the past as expected, but at redshifts z ~ 1, w drifts from its GR value and 
presents phantom behaviour, w < —1. We note that the size of this 'drift' depends almost 
entirely on our choice of b; the larger we take 6, the smaller R/e is and the more suppressed 
the tanh(i?/e — b) and sech^(i?/e — b) terms become. This is clearly shown in fig. 12; as we 
increase b, the late time drift from w = —1 becomes increasingly suppressed; \Aw\ ~ 0.06 
for b = 1.2 and \ Aw\ ~ 2 x 10~^ for 6 = 4. 

Finally, we note that w crosses the phantom boundary and satisfies w > —1 for z < 1 
and moderate values of b. Thus, our model naturally exhibits phantom behaviour at recent 
redshifts during the matter dominated stage. Deviations of w{z) from —1 are small, less 
than several percent, hence there is a good agreement with present observational upper 
bounds on \wbe + 1|i see e.g. [6]. 



w 




Figure 12: The 'dark energy' equation of state parameter w as a function of redshift for the gR^- 
AB model, with b = 1.2 (top) and 6 = 4 (bottom). We note that w oscillates around w = —I at 
earlier times, as can clearly be seen in the 6 = 4 case, however for redshifts z ~ the oscillations are 
small but the 'drift' terms in p and Pp become important. We see the deviation at z ^ depends 
on 6; as 6 increases, deviations from GR become increasingly suppressed. 
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5. Conclusions and discussion 

To summarize, in this paper we have shown for the first time that it is possible to construct 
at least one self-consistent model of present DE in the scope of F{R) gravity which sat- 
isfies all five viability conditions presented in the Introduction, is free of new singularities 
and does not destroy any of the previous successes of cosmology. To achieve this aim, it 
has been necessary to extend the range of R over which previous DE F{R) models were 
defined to both large positive and negative values, and change the behaviour of the models 
correspondingly; see Eq. (4.2) for the improved gR^ -AB model. 

Furthermore, since the large- i? behaviour (3.1) of the F{R) function, which is needed 
to avoid new singularities and to solve two other problems of previous models, appears to be 
just the same as needed for scalaron-driven inflation in F{R) gravity, wc have shown that 
the model Eq. (4.2) can describe inflation (primordial DE), the present acceleration of the 
Universe (present DE) and the intermediate epochs of radiation and matter domination for 
the unique choice of its parameter M, determined by the observed power of scalar (density) 
perturbations. Unexpectedly, we have found that the low-curvature modification of F{R) 
from its GR value, which is needed to describe present DE, strongly affects processes at very 
high values of curvature, specifically during reheating after inflation, through its non-zero 
^-factor (4.1). As a result, in contrast to pure inflationary models of F{R) gravity which 
have g = 0, scalaron oscillations after the end of inflation become strongly non-linear 
and the Universe evolution passes through the sequence of interchanging periods with 
a{t) oc \/t and a const, with the number of time e-folds Alnt equal to 21n(l/(l — 2^)) 
and ln(l/(l — 23")) for them correspondingly. We find that on average, the Universe expands 
as a{t) oc i^/^ during this period; this is the most significant mathematical result of the 
paper. 

Creation of particles and antiparticles of usual matter and final reheating axe achieved 
by taking into account the process of gravitational particle creation by these non-linear 
oscillations of R. This mainly occurs at the end of inflation, so reheating is less efficient 
than in the purely inflationary R + R'^/QM'^ model, although still viable. Due to the 
different average law of expansion after the end of scalaron-driven inflation, predictions for 
parameters of primordial spectra of scalar (density) perturbations and gravitational waves 
generated during inflation in this combined model of primordial and present DE are slightly 
different to those for the inflationary model only. The difference is due to the change in the 
number of scale factor e-folds N used in the corresponding formulas: from N ~ (50 — 55) 
to = 70. This specific prediction of the combined model is observable, however, it is 
degenerate with a possible slow variation of the R^ behaviour of F[R) at large R in the 
model (4.2), e.g. like in Eq. (4.6). 

In this paper we have not addressed the issue of neutron star stability in F{R) gravity 
raised in [47] and further considered in [59, 57], since the most recent results in [60, 61] 
suggest that there is no problem and that, for our choice of the parameter M in the 
improved model (4.2), non-GR corrections are very small if the trace T of the matter energy- 
momentum tensor remains non-positive inside neutron stars. However, due to a non-zero 
^(-factor, this problem becomes highly non-trivial and requires special consideration if T is 
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permitted to become positive at large matter densities. 

The most critical observational prediction of DE models in F(R) gravity remains the 
anomalous growth of density perturbations in the matter component at recent redshifts 
[44, 71, 72, 73, 74] which results, in particular, in a mismatch between parameters such as 
as and Ug determined from CMB temperature anisotropy and galaxy clustering separately, 
assuming GR (as well as from the cluster abundance at different z). The absence of such 
an effect at the level of ~ 5% in the HSS model makes its background evolution practically 
indistinguishable from that in the standard ACDM. However, this effect is suppressed in 
the AB model as compared to the HSS one, so it is more difficult to falsify the former 
model. Future observational data will determine the fate of this whole class of models. 
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